Static and Dynamic Critical Behavior of a Symmetrical Binary 

Fluid: A Computer Simulation 

Subir K. Das/ Jiirgen Horbach,^ Kurt Binder,^ Michael E. Fisher/ and Jan V. Sengers^ 

^ Institute for Physical Science and Technology, 
University of Maryland, College Park, MD 20742, USA 
^ Institut fiir Physik, Johannes Gutenberg Universitdt Mainz, 
Staudinger Weg 7, 55099 Mainz, Germany 
(Dated: February 6, 2008) 

Abstract 

A symmetrical binary, A+B Lennard-Jones mixture is studied by a combination of semi- 
grandcanonical Monte Carlo (SGMC) and Molecular Dynamics (MD) methods near a liquid-liquid 
critical temperature Tc- Choosing equal chemical potentials for the two species, the SGMC switches 
identities (A — > B ^ A) to generate well-equilibrated configurations of the system on the coex- 
istence curve for T < Tc and at the critical concentration, Xc = 1/2, for T > Tc- A finite-size 
scaling analysis of the concentration susceptibility above Tc and of the order parameter below is 
performed, varying the number of particles from N = 400 to 12800. The data are fully compatible 
with the expected critical exponents of the three-dimensional Ising universality class. 

The equilibrium configurations from the SGMC runs are used as initial states for microcanonical 
MD runs, from which transport coefficients are extracted. Self-diffusion coefficients are obtained 
from the Einstein relation, while the interdiffusion coefficient and the shear viscosity are estimated 
from Green-Kubo expressions. As expected, the self-diffusion constant does not display a detectable 
critical anomaly. With appropriate finite-size scaling analysis, we show that the simulation data for 
the shear viscosity and the mutual diffusion constant are quite consistent both with the theoretically 
predicted behavior, including the critical exponents and amplitudes, and with the most accurate 
experimental evidence. 

PACS numbers: 47.27.ek, 64.60.Fr, 66.10.Cb, 64.60.Ht 
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I. INTRODUCTION 



Recently, there has been a renewed interest in the critical behavior of simple and com- 
plex fluids, both with respect to liquid-gas transitions and demixing transitions in binary 
fluids . There has been remarkable progress both in the fuller theoretical under- 
standing of the asymptotic critical region and the types of correction to scaling that occur 
in these systems, -"^ and of the nature of the crossover towards classical, van der Waals-like 
behavior further away from the critical point. ^ Also, the values of the critical exponents and 
other universal properties (such as critical amplitude ratios) are now known with high preci- 
sion, from a variety of techniques (renormalization group,^ Monte Carlo simulations,^^ and 
high-temperature series expansions^). Particularly relevant in the present context are ad- 
vances in the flnite-size scaling analysis of computer simulations of fluids,'i'^'2iiii which allow 
one to study both universal and non-universal critical properties of various off-lattice models 
of fluids with an accuracy that is competitive with the work on Ising lattice models.^^ 

With respect to critical dynamics in fluids, the situation is much less satisfactory even 
though precise experimental data were presented a long time agc^^ and theoretical analy- 
ses invoking either approximations, such as mode coupling theor y -^^i-^^i-^'^i-^^i-^^i-^''' or low-order 
renormalization group expansions in e = 4 — d, where d is the dimensionality,— i^^iiSiSi do ex- 
ist. One should note that dynamic universality classes encompass fewer systems than static 
ones j'-~i-- while uniaxial ferromagnets, binary alloys, liquid-gas criticality and demixing in 
binary fluids all belong to the same universality class as far as their static critical behavior 
is concerned, these systems belong to more than one dynamic universality class. Thus, there 
is a clear need for more theoretical analyses of critical dynamics. 

Particularly scarce are simulations of the critical dynamics of fluids: of recent works, 
we are aware only of a nonequilibrium Molecular Dynamics (NEMD) calculation for a two- 
dimensional fluid in which heat conduction near the critical point was studied^^^ and of a 
similar investigation by Chen et ali^ of a three-dimensional Lennard- Jones single component 
fluid. Beyond that Jagannathan and Yethiraj (JY)^ used a three-dimensional Widom- 
Rowlinson model^ to study the inter-diffusional critical dynamics in a binary-fluid. However, 
the conclusions of this latter work have been seriously challenged.— 

In the present paper, we move to flll this gap by presenting a comprehensive simulation 
study of critical dynamics in fluids by studying a symmetric binary fluid Lennard- Jones 
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mixture. In previous work2LS2i22 we have shown that the coexistence curve, concentration 
susceptibihty, interfacial tension between coexisting hquid phases, pair correlation func- 
tions and static and dynamic structure factors for this model can be reliably estimated 
via a combination of semi-grandcanonical Monte Carlo methods (SGMC) ^'^i^^i^^i^'^i^^i^^i^^i^^ 
and microcanonical Molecular Dynamics (MD) methods .^^i^^i^'^'i^^ Transport coefficients such 
as self-diffusion and interdiffusion coefficients,— shear viscosity^ and bulk viscosity^ were 
estimated away from the critical region. Here, however, we expressly address the critical 
behavior of the model and compare with theoretical expectations. Moreover, our study 
strongly supports the challenge^^ to the earlier study by JY— of the somewhat similar but 
less realistic Widom-Rowlinson model. 

In the balance of this article. Section II presents the model and briefly reviews our sim- 
ulation methods. Section III discusses the static critical properties that are extracted from 
the "raw data" by a finite-size scaling analysis. ^^^iiML^ Sec. IV then presents our results 
on the selfdiffusion coefficients and the shear viscosity: we discuss them in the light of the 
theoretical predictions .^'^' ' ^^^^^' ' ^^^^^' ' ^'^^^^' ' ^^' ' '^^' ' '^^ The interdiffusional coefficient, which vanishes 
fairly strongly, requires more detailed, specifically, finite-size scaling considerations, etc. 
which are presented in Sec. V. The article concludes in Sec. VI with a brief summary and 
discussion. 



II. MODEL AND SIMULATION TECHNIQUES 

Following Refs. 27 — 29 we consider a binary fluid of point particles with pairwise inter- 
actions in a cubical box of flnite volume V = subject to periodic boundary conditions. 
Starting from a full Lennard- Jones potential 

*Lj(rij) = 4e^/3[(a^/3/rij)^2 - [a^^/rijf] , (1) 

we construct a truncated potential that is strictly zero for Vij > Tc as follows,— 

u{rij) = $Lj(rij) - ^Lj{rc) - {rij - rc)^^^^^\r,,=r,, for Vij < r^. (2) 

This form ensures that both the potential and the force are continuous at r = r^. In the 
previous work,i^>2Si2£ the last term on the right-hand-side of Eq. Q was not included so that 
the force at Vij = exhibited a jump, while only the potential was continuous. This is not 
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desirable when considering dynamic behavior, because, in a microcanonical simulation, this 
results in a drift of the total energy. 

The parameters in Eqs. dH) and © were chosen as 

0"AA = O-RB = CTAB = 0- , (3) 

and, hence, we adopt a as our unit of length. The cutoff is chosen as = 2.5cr. As 
previously,— we take the total particle number = Aa + A'b of the binary A+B mixture 
and the volume such that the density p* = pa^ = Na^/V = 1. This choice is convenient, 
since the system is then in its liquid (rather than its vapor) phase and crystallization is not 
yet a problem at the temperatures of interest. Finally, the reduced temperature T* and 
energy parameters Eap are chosen as2L22i2S 

T* = kBT/e and ^aa = £^bb = 2£:ab = (4) 

The system is equilibrated as follows. First, a Monte Carlo (MC) run is performed 
in the canonical ensemble {N\=N-b,V,T), starting out from particles at random posi- 
tions in the simulation box. The MC moves used are random displacements of randomly 
chosen single particles (selecting the trial value of each new cartesian coordinate in the 
range [— cr/20, +a/20] about its old value) and applying the standard Metropolis acceptance 
criterion . ^^1^^ These initial runs were carried out for 10^ Monte Carlo steps (MCS) per par- 
ticle, for systems with particle numbers from A^ = 400 to A^ = 12800. Then equilibration 
is continued, using the semi-grandcanonical Monte Carlo (SGMC) method . ^'''i^^i^^i^^i^^i^^i^^i^^ 
After 10 displacement steps per particle A^/10 particles are randomly chosen in succession 
and an attempted identity switch is made, A B or B — A, where both the energy change, 
AE, and the chemical potential difference, A/i, between A and B particles enters the Boltz- 
mann factor. However, we restrict attention here to the special case Ap = 0, which, for the 
symmetric mixture considered, means that for T > T^. we simulate states with an average 
concentration (a;A) = {xb) = 1/2 (with Xa = Na/N), while for T < Tc we simulate states 
along either the A-rich or the B-rich branch of the coexistence curve. 

In the semi-grandcanonical ensemble, the concentration xa is a fluctuating variable, and 
hence the probability distribution P{xa) can be recorded: this is particularly useful in the 
context of a finite-size scaling analysis.— >^2iiL^ In addition, use of the semi-grandcanoncial 
ensemble for the study of static critical properties in a binary fluid mixture is preferable 



since critical slowing down is somewhat less severe than in the canonical ensemble. Critical 
slowing down limits the accuracy that can be attained, since the statistical error scales like 
1 / y^, where n is the number of statistically independent configurations.^— The statistically 
independent configurations used in the computation of averages must be separated from 
each other along the stochastic (MC) or deterministic (MD) trajectory of the simulation by 
a time interval which is of order of the longest relaxation time in the system.— In a finite 
system at the demixing critical point of a fiuid binary mixture, this slowest relaxation time 
scales with box dimension L as 

(5) 

where 2; is a universal dynamic critical exponent, which depends on the dynamic universal- 
ity class For the SGMC algorithm, the order parameter (concentration difference 
between A and B) is not a conserved variable, while the average density p = N/V is con- 
served. As a consequence, this model belongs to "class C " in the Hohenberg-Halperin 
classification,^ and hence the dynamic exponent is roughly z = If we performed MC 

simulations in the NaN^VT ensemble, the order parameter would be a conserved variable, 
in addition to the density ("class D"-^), and then the dynamic critical exponent is expected 
to be significantly larger, z = 4 — t], where rj (^ 0.03)^»^ describes the spatial decay of 
correlations at criticality. If one performs MD runs in a microcanonical A^aA^b^-^' ensemble 
(with Aa = A^B and with an energy chosen so that the system is precisely at the critical 
point), the dynamic exponent is predicted to be 2; ~ 3 _i6ii7^i8 fact, Jagannathan and 
Yethiraj^ used Eq. dSJ in order to explore the critical dynamics of their Widom-Rowlinson 
model.-^ Comparing the values of z for the three algorithms discussed above, it is clear why 
the SGMC algorithm has an advantage, as far as static critical properties are concerned. 

For the study of dynamic critical properties, multiple independent initial configurations 
were prepared,— from SGMC runs with 5 x 10^ MCS (but excluding states from the first 
10^ MCS). Then a further thermalization for 2 x 10^ MD steps was carried out in the NVT 
ensemble using the Andersen thermostat .■^^'^^^^^'^ For all MD runs, we always chose the 
masses of the particles equal to each other, ttia = ttlb = m, and applied the standard 
velocity Verlet algorith m^^i^^i^'^'i^^ with a time step 6t* = O.Ol/v/48 where t* = t/to with 
scale factor 

to = {ma'/ey/'. (6) 
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The final production runs in the microcanonical ensemble (where the Andersen thermostat 
is switched off) used from about 10^ MD steps for temperatures T* = 1.5 and higher, but 
up to 2.8x10^ MD steps for T close to (where T* = 1.4230 ± 0.0005, see below). To 
avoid confusion, we note that the different value of in the previous work,"2L2Si2S namely, 
T* = 1.638 ± 0.005, arose from the different choice of potential (the last term on the right 
hand side of Eq. Q being absent in Refs. 27 — 29). 

III. STATIC CRITICAL PROPERTIES 

Using the SGMC algorithm we record the fluctuating number of A particles A^a (recall 
that A^B = N — Na with the total particle number A^ fixed) and generate histograms to esti- 
mate the probability distribution P{xa) of the relative concentration xa = Na/N. Typical 
"raw data" for P{xa) are shown in Fig. ^ The symmetry relation that holds for A/i = 0, 
namely, 

P(xa) =P(1-xa), (7) 

has been incorporated in the data. This has been done because below Tc, where P{xa) has 
a pronounced double-peak structure corresponding to the two sides x™'^''''^'* of the 

two-phase coexistence curve, transitions from one side to the other occur very infrequently 
(or, at low temperatures, not at all). In Fig. ^a) we present probability distributions for 
several temperatures below Tc while Fig.^b) shows the distributions for temperatures above 
Tc- From P{xa), we define the truncated moments (a;^) of the concentration distribution as 
follows^- 

(a;A) = 2 / XAP{xA)dxA , (8) 

Jl/2 

(xi) = 2 / xiP{xA)dxA . (9) 

Jl/2 

The two branches of the coexistence curve can then be estimated for large A^ via 

xr(^)^(l-(a:A)), xr^'^c^ixA), (10) 

while the "concentration susceptibility" x and its dimensionless form, can be estimated 
from 

kBTx = T*x* = Ni{xl)-l/4), T>Tc. (11) 
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Another useful quantity is the fourth-order cumulant Ul, defined bj^^^^Sdi 

U,{T) = ((XA - l/2)^)/[((xA - 1/2)^)2]. (12) 

Note that for a finite system, {{xa) — 1/2) remains nonzero even for T > Tc [as seen 
in Fig. ^b)]. Furthermore (xa) is a smooth function of temperature for finite L while x 
likewise remains finite at Tc. Due to these effects, a finite-size scaling analysis of the SGMC 
data for these quantities is clearly required, as is well known.— 

There are different strategies used in the literature to estimate the location of the critical 
temperature Tc from such simulation results . The simplest method, most 
often used for fluids, records x™'^^*'^^ , x^*'^*"^'' for several choices of L and checks for a regime 
of temperatures below Tc where the results are independent of L within statistical errors. In 
this regime one chooses several temperatures, as close to Tc as possible, and fits to a power 
law 
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where the critical amplitude B and Tc are adjustable parameters, while the critical exponent 
(3 is fixed at its theoretical value for the universality class of the three-dimensional Ising 
model, P ^ 0.325.-2'^ 

In Fig. El we show the two-phase coexistence curve for N = 6400. (Recall that for our 
choice we have density p* = 1, so the system size is L = N^^^a.) The dashed line in Fig. |21 
is a guide to the eye for the numerical data (open circles). The continuous line is a fit to 
Eq. (fT^ using the range from 0.2 < xa < 0.5 (but excluding the two points closest to 
Tc). From this we obtain T* = 1.423 ± 0.002 and B = 1.53 ± 0.05. This fit is good over a 
relatively wide range of temperature but, in reality, the range over which Eq. ()13|1 should 
be valid is significantly smaller owing to various corrections to scaling^*^^ which have been 
neglected. Hence, systematic uncontrollable errors easily arise: the true value of Tc could 
well be somewhat lower with B larger. Alternative methods are thus needed to obtain 
more reliable estimates of Tc with improved error bounds. Indeed one clear drawback of 
the previous study on critical dynamics by JY~ is that the accuracy with which the critical 
point could be located was relatively poor. 

A method used more recently in Monte Carlo studies of critical phenomen a?'^''^^'^^ is data 
collapsing based on the finite-size scaling hypothesis.—"^ Let us consider the concentration 
susceptibility x i^i the vicinity of the critical point. At the critical concentration, xa = xb = 



1/2, we have 

X*(T)^roe-^ with e = {T-n)/T,, (14) 

where Fq and 7 are the critical amphtude and critical exponent, respectively. From here on 
we shall use the symbol e for the reduced temperature deviation to avoid confusion with t 
which arises naturally as a symbol for time in the context of MD simulations. For a system 
of finite size L, the basic scaling ansatz may be written as 

xliT) ^ roZ{y)e-\ (15) 

where y = L/C,, in which ^ is the correlation length, while Z{y) is the appropriate finite-size 
scaling function. The correlation length diverges at criticality as 

e ^ eoe^^ (16) 

with amplitude ^0 and exponent u. In the limit y ^ 00 (so that L — 00 at fixed e > 0) 
the scaling function Z{y) must approach unity so that Eq. (|14p is recovered. For static 
quantities (such as x)^ iii short-range systems with periodic boundary conditions, Z{y) 
generally approaches unity exponentially fast. So, one may expect the behavior 

Z{y) = l + Z^y^e-"y + ..., for y ^ 00, (17) 

where the values of the exponent ip and the integer n = 1, 2, 3, ... depend upon the details 
of the system in question. On the other hand, for finite L, in the limit y ^ (so e — at 
fixed L < 00), the susceptibility Xl{T) is finite and its variation with T must be smooth 
and analytic. Thus one should have 

Z{y)=f'/''[Zo + Z,y'/'' + Z,y'/'' + ...] as y ^ 0. (18) 

An effective procedure is then to study the computed quantity x|^(T)e'' as a function 
of y for different system sizes by using Tc as an adjustable parameter to optimize the data 
collapse. The plot should then approach roZ{y). From the asymptotic behavior of ToZly) as 
y becomes large one can then estimate the critical amplitude. In Fig. |21we plot x*l{T)^'' vs. 
y/{y + yo) noting that the abscissa variable approaches zero when y but tends to unity 
when ?/ — » 00 for yo > 0. For convenience we have chosen yo = 7 which is comparable with 
the value of L/^ for the largest system size. (This point will be discussed further.) Of course, 
if one wishes to estimate all three quantities T^, 7 and u from such a procedure, one again 
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fights hard to control systematic errors since the uncertainities in the estimates for these 
quantities are inevitably highly correlated. Accordingly we have fixed the critical exponents 
at their universal values, accurately known from other studies,— 7 ~ 1.239, v ~ 0.629, 
since there is no reason to doubt that these exponents describe the static critical behavior 
of the present model. 

In Fig El we demonstrate data collapse for four trial values of Tc. It is clear that the 
collapse is inferior for the values = 1.425 and 1.419 compared to the choices 1.423 and 
1.421: the collapse looks quite acceptable in these latter cases. Thus, we conclude that our 
previous estimation of Tc is consistent with this analysis. 

An unbiased method to evaluate utilizes the fourth-order cumulant defined in 
Eq.( IT^ :^°i'^^ plotting Ul{T) vs. T for various sizes L one finds Tc from a common inter- 
section point of these curves once corrections to finite-size scaling become negligible. Quite 
generally one has — 1/3 in any one-phase region while Ul 1 on the coexistence-curve 
diameter in the two-phase region. Furthermore, for the three-dimensional Ising universality 
class Ul{Tc) takes the value 0.6236..2*a 

In Fig. ^a) we present Ul{T) for several system sizes (as indicated in the figure) over 
a rather wide range of temperature. The horizontal dashed line indicates the value of this 
quantity at the critical point for the Ising universality class. This plot clearly confirms that 
our model belongs to the three-dimensional Ising universality class. While on a coarse scale 
the expected intersection is nicely seen, the enlarged view of the data, in FigEJ^b), reveals 
some scatter, which is mainly due to the statistical errors of the simulation data [as can be 
seen, by comparing with plots in Refs. 8, 9(c) and 45]. In light of these statistical errors, 
further analyses are clearly not warranted. Thus, for example, the method of Wilding^ based 
on the use of the full distribution P{xa) at criticality, rather than using only the second 
and fourth moment, is not tried here: but see also the critique in Ref. 46. In Fig. |^b), the 
smooth lines are fits to the hyperbolic tangent function. For the system sizes shown, these 
fits all intersect one another close to T* = 1.423 and at the Ising critical valued 0.6236. Since 
this method appears to be the most reliable currently available in the literature, we adopt 

T* = 1.4230 ± 0.0005 (19) 

for the subsequent analysis of our simulation data. 

Of course, it is also of interest to extract the pair correlation functions gAA{r), gAB{r) 
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and gBsir) and their Fourier transforms, S'aa('?), Sab{(i) and 5'bb('?), as described for a 
closely related model (outside the critical region) in earlier work.— ^ Of particular interest 
are the combinations that single out number-density fluctuations Snnil) and concentration 
fluctuations Scdq), defined via^ 

Snniq) = ^AA(g) + 2^AB(g) + ^BB(g) , (20) 

S'cc(g) = (1 - XA^SAAiq) + xj^SBBiq) - 2a;A(l - XA)SABiq) ■ (21) 

Fig- El^) shows that Snn{q) exhibits the normal oscillatory behavior of the structure factor 
of a dense liquid: the approach to criticality has little discernible effect. By contrast, 
Scc{q) varies weakly at large q, corresponding to small spatial length scales, while at small 
q a strong increase is observed. Of course, this is the expected Ornstein-Zernike behavior 
reflecting the "critical opalescence" due to concentration fluctuations: see Fig. Efb). The 
inset here displays an Ornstein-Zernike plot, based on 

SUq) = kBTx/[l + q'e + ■■■], (22) 

from which our estimates of the correlation length ^ have been extracted. Note that, in the 
fitting process, we have used the value of kBTx from our SGMC simulations. 

In Fig. El we plot the susceptibility and correlation length x*{T) and C,{T), as functions 
of e for = 6400 and fit the data with the respective asymptotic forms (jl4j) and For 
the fitting we again adopt the Ising critical exponent values, so that the amplitudes are the 
only adjustable parameters. The quality of the fits suggests that the finite-size effects are 
negligible in this temperature range (where as one sees from Fig. El ?/ > 4 so that L > 4^(T)). 
The amplitudes are found to be 

^o/a = 0.395 ± 0.025, Fq = 0.076 ± 0.006. (23) 

In Fig. El^c) we plot x vs. ^. The continuous line is a power-law fit with the exponent 
7/z/ ^ 1.970. 

We conclude this section by noting that no unexpected features have been uncovered. 
The static properties comply fully with the anticipated critical behavior characterizing three- 
dimensional Ising-type systems. Furthermore, the corrections to scaling seem to be quite 
small in the temperature range covered by our simulations, so that one can observe the 
asymptotic power laws even relatively far from the critical temperature. 
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IV. SELF-DIFFUSION COEFFICIENT AND SHEAR VISCOSITY NEAR CRIT- 
ICALITY 

The transport coefficients which are most readily found from simulations are the self- 
diffusion coefficients and which can be extracted from the Einstein relations for the 
mean square displacements of tagged particles,-^ namely, 



and likewise for gB{t), where it is understood that the average (. . .) includes an average over 
all particles of type A or B, respectively. The self-diffusion coefficients then follow from 



and similarly for Db- In the region above Tc and at critical concentration xa = = 1/2, to 
which we will restrict attention, the symmetry of the model requires gA{t) = gsit) and 
Da = Db = D. To extract D and study its temperature dependence we have hence 
averaged over the mean square displacements of all particles: see Fig. |7| Note that in 
the initial, ballistic regime gA{t) varies quadratically with t before crossing over to linear 
diffusive behavior from which Da is estimated. As expected, the temperature dependence 
of D is rather weak and, indeed, close to linear over this fairly narrow temperature range; 
moreover, D appears to remain nonzero at T = T^, with a value close to D* = 0.048. Indeed, 
there is no sign of any critical anomaly, consistent with the previous work of JY.— Similarly, 
a study of self-diffusion near the vapor-liquid critical point of a lattice gas model^ did not 
detect any significant critical anomaly. (Note, however, that this latter model belongs to 
class B in the Hohenberg-Halperin classification.i^) Nevertheless, a weak anomalous decrease 
of the self-diffusion coefficient at the critical density has been seen in simulations of simple 
fluids near the vapor-liquid critical point, ^'^ but has not yet been conffimed experimentally.— 
Next we consider the reduced shear viscosity rj* which we calculate from the Green-Kubo 
formula!^ 



gAit) = ([f.,A(o)-n,A(t)]'), 



(24) 



D*A = ito/a^)DA = (to/a') \im[gAit)/6t] 



(25) 




(26) 



where the pressure tensor Cxyit) is given by 




(27) 



i=l 
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in which Vi is the velocity of particle i while F is the force acting between particles 
Theoryi^»i^^ suggests that at the critical point rj* should diverge as 

V* ^ Vo^'"'"'' - C\ (28) 

where rjo and > are the appropriate critical amplitude and exponent. Renormalization- 
group theorj*iS»i^ predicts a;^ ^ 0.065 while the theory of Ferrell and Bhattacharjee^ yields 

^ 0.068. These values are consistent with experiments^!^ which yield between 0.064 
and 0.070. The most recent theoretical estimate is^ 0.0679 ± 0.0007 and the most recent 
experimental value is^ 0.0690 ± 0.0006. 

Fig.lHldisplays a log- log plot of rj* vs. e. As normal in MD simulations, accurate estimation 
of the shear viscosity is difficult and the large error bars shown in Fig. |H1 prevent us from 
making definitive statements about the critical singularity. But the slow increase of rj* as 
T ^ Tc is, in fact, compatible with the expected power law divergence, as Fig. |H1 shows, 
since the fitted line has a slope corresponding to u = 0.629 and Xr^ = 0.068 in Eq. (f28|l . 
Although this fit is consistent with the theoretical prediction, estimating i/x^ from the data 
itself is clearly of little value. However, the amplitude, for which we obtain rjo = 3.87 ± 0.3, 
will prove to be significant. 

At this point it is of interest to check the validity of the Stokes-Einstein relation, which 
relates the self-diffusion constant of a diffusing spherical particle of diameter d to the shear 
viscosity of the fluid ri{T). For a particle moving in a fluid of like particles, the Stokes- 
Einstein diameter d can be written^! 

d = aT*/2nri*D*, (29) 

which corresponds to the assumption of slip boundary conditions on the surface of the 
diffusing sphere. For stick boundary conditions, a factor 3 replaces the factor 2 in Eq. ()29j) . 
In Fig. ini we present a plot of d in the interval T* = 1.45 to 1.55. The data suggest that 
Eq. fj29p is still a valid approximation despite the strong concentration fluctuations close 
to Tc- However, we do not expect the relation to remain valid closer still to since 1]{T) 
diverges, albeit slowly, while D remains flnite. 
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V. INTERDIFFUSION: FINITE-SIZE SCALING 



Finally we consider the interdiffusion or mutual diffusion coefficient, which is expected 
to vanish at the critical point. Following the previous work^L^SiSS we use the Green-Kubo 
formulaF^ which we express as 

2 2 

DMT) = ^DUT) = f hm D1b(T; t), (30) 
where, introducing the appropriate reduced Onsager coefficient, 

C{T) = Urn C{T;t), (31) 

t— >oo 

we have the relation 

DlB(T;t) = £(T;t)/x*(T), (32) 
in which, for numerical purposes, we will use our fits to Eq. (fT^ for x*, while 

£(T; t) = (to/Na'n f dt'{J^^{Q)J^^it')), (33) 
in which the current vector J^^ is defined by 

jAB (^) = ^ ^. ^ (^) _ Jr^Y,V^,B{t) ^ (34) 
i=l i=l 

where Vi^aif) denotes the velocity of particle i of type a at time t. 

Since the correlation function (</^^(0) J^^(t)) is rather noisy at large times comparable 
to the total observation time, tobs, of an MD run, it is difficult to attain a precision for 
Z^AB as high as for the self-diffusion constant D. This problem is evident in Fig. where 
short-time maxima (seen for t* ~ 0.1) are followed by shallow minima (at t* ~ 0.5 — 0.6) 
that appear before the expected plateau starts at about t* ~ 4. For times t* > 40, the curves 
become progressively more noisy, and, clearly, at most temperatures the data for t* > 100 
may be discarded owing to deficient statistics. But, in any case, the facts that the general 
shape of the curves is similar for all temperatures studied, and, in particular, that the time 
needed for DAB(t) to reach a plateau value is almost independent of T, suggest that the 
dominant contribution to the temperature dependence of -Dab arises from x* in Eq. (j32j) . 

Fig. 111! presents plots of D\-q{T) versus T for systems of = 6400 particles. One sees 
that -Dab has the apparent power-law behavior 

D'ab ~ C'"^*' ~ e^'^""" ^ e\ with Xesu ~ 1, so Xes ^ 1-6. (35) 
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This result is in strong disagreement with the theoretical prediction for the interdiffusional 
critical exponent,— namely, 

= 1 + ~ 1.0679, (36) 

that should be accessible asymptotically when T — > T^. Indeed, the value of Xes found from 
Fig. ^2 is even larger than the value xd = 1-26 ± 0.08 quoted by JY— for their binary fluid 
model. However, it would be quite erroneous to conclude from our data and Fig. ^2 that 
the simulations indicate a serious failure of the theory. Even though we have found that the 
finite-size effects in the equilibrium static properties are small for the temperature range and 
system sizes explored [where L > 4^(T)], one must be prepared to encounter much stronger 
finite-size corrections in transport properties near criticality. 

In addition experiments have shown that the Onsager coefficient near a liquid-liquid 
demixing critical point or its equivalent, the thermal conductivity near a vapor-liquid critical 
point, may have a significant noncritical background contribution arising from short-range 
fluctuations.-^^ 

Thus it is crucial to analyze dynamical simulations by making proper allowance for the 
finite-size behavior and also for possibly significant 'background' contributions to the quanti- 
ties computed. To this end we consider, as above, only the critical isopleth xa = = 1/2 
and will focus on the finite-size Onsager coefficient Cl{T) as defined via Eq. (jHS))- The 
prime reason for analyzing the Onsager coefficient rather than the interdiffusional constant 
Da_b{T) (which clearly follows by dividing by x) is that it represents most directly the basic 
fluctuation sum/integral analogous to expressions like Eq. (jlip defining or the standard 
fluctuation sums for the specific heat, etc.;— experience teaches that such properties display 
the simplest, albeit not "simple", singularity structure and finite-size behavior. 

The variation of Cl{T) with temperature (on a linear scale) is presented in Fig. IT^ for 
the largest computationally feasible system-size of = 6400 particles and, thus, of box size 
L ~ 18.6(T. Note, first, that although C(T) rises sharply close to Tc, there seems to be a 
relatively large background contribution. If one chose to ignore this and merely examined 
a direct log-log plot extending over only one decade in e, the resulting effective exponent 
would be of little theoretical significance. As shown by Sengers and coworker s,^^'^^*^ the 
Onsager coefficient close to criticality may be written (in the thermodynamic limit, L —>■ oo) 
as 

C{T)=C,{T) + AC{T), (37) 
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where Cb{T) is a slowly varying background term which arises from fluctuations at small 
length scales, of order a, while AC{T) represents the "critical enhancement" induced by 
long-range fluctuations on the scale of the diverging correlation length ^(T) ~ e~'^ . The 
singular piece is predicted to diverge as 

A£{T) « QTye"' (38) 

with an amplitude Q and an exponent 

= x\u ~ 0.567 with xa = 1 — r/ — , (39) 

where is defined via Eq. ()28|) and expresses the weak divergence of the viscosity ^^(T), while 
Tj ~ 0.03 is the standard critical exponenlj^ [that enters the scaling relation 7 = (2 — ini)v\. 
Furthermore, we may invoke the "extended Stokes-Einstein relation" ^^1^° 

AL>AB(r) ^ RDkBT/67cr]{T)aT), (40) 

which embodies the relation xd = 1 + [see and Ref. 11]. This is expected to describe 
the singular part, ADab, of the mutual diffusion coefficient Dab{T) and so leads to the 
explicit expression 

Q = RdTo (T/67Tr]o^o (41) 

for the amplitude in (jHHj) . Here Rd is a universal constant of order unity while To, tjq, and 
^0 are the critical amplitudes for x*, 77*, and ^ defined via Eqs. (IHj), (jlHl) and (fTHj). 

Two theoretical methods have been developed to calculate the universal dynamic 
amplitude ratio Rn, namely, mode-coupling theory of critical dynamics and dynamic 
renormalization-group theory.^^ In first approximation, mode-coupling theory yields^^ Rd = 
1.00; but when memory and nonlocal effects are included one obtains the improved estimate^ 
Rd = 1.03. The early theoretical values obtained from renormalization-group theory have 
varied from 0.8 to 1.2 due to various approximations, as reviewed by Folk and Moser.— 
The calculation of Folk and Moser with the fewest approximations has yielded Rd = 1.063. 
Experiments give values consistent with the mode-coupling prediction . Here we follow 
Luettmer-Strathmann et al}^ and adopt the estimate Rd = 1.05 as a compromise between 
the predictions of mode-coupling theory and the renormalization-group calculations. Using 
the estimates reported above for the amplitudes in (PT|) then yields 

Q = (2.8 ± 0.4) X 10"^ (42) 
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numerical prediction for the present model. 

Our aim now is to discover if this theoretical analysis and the value ()42j) for Q are 
consistent with the evidence available from our simulations. Because of the computational 
demands imposed by the collective transport properties we have obtained results over a 
substantial temperature range only for the (A^ = 6400)-particle system; however, for T* = 
1.48 ~ 1.04 T* we also computed Cl{T) for = 400, 1600, and 3200. 

To analyze these data we write the finite-size scaling ansatz as^MSiS 



where y = L/C, while W{y) is the finite-size scaling function. As already discussed in the 
context of static critical phenomena, one requires W{y) —>■ 1 when y ^ oo so as to reproduce 
the correct behavior (jHHj) in the thermodynamic limit. In this case, however, it is not clear 
how rapidly W{y) should approach unity. Indeed, since transport properties are calculated 
from time correlation functions of currents, they reflect the nonequilibrium behavior of the 
system. Although the exponentially rapid approach that applies in the static case [see 
Eq. p7|) ] might still be realized here, the well known, noncritical long-time tails in the 
correlation functions, etc., suggest that a slower, power-law approach cannot be excluded. 

On the other hand, for finite L all properties remain bounded through criticality so that 
in the limit y —* one should have 



as in (fTHj) . where from (jHIJj) we have xx — 0.90. 

Of course, we do not know the value of the background Cb(T) in but since it is 

slowly varying, we may reasonably replace it by a constant effective value Cl^^ . Then, by 
treating Cl^^ as an adjustable parameter and examining the simulation data for >Vl(T) = 
A£i(T)e'^^/T* as T and L vary with ux set to its Ising value, we may seek an optimal data 
collapse onto the scaling form QW[L/C(T)]. Note that if this is achieved, the value Q should 
emerge when y = L/C, becomes large.— 

Fig. 1121 presents separated plots of Wl(T) vs. y/{y + y^), with, as in Fig|21 yo = 7, for 
four assignments of C^^^ . Note that the filled symbols represent the data at T* = 1.48 for 
system sizes L/cr ~ 7.37, 11.70, 14.74, and 18.57; their reasonably good collapse onto the 
remaining data (all for L/a 18.57) and their approach towards for small y serve to 



ACl{T) = C{T) - C,{T) ^ Q T*W{y)/e''\ 



(43) 



W{y)^Woy'''[l + Wiy^l^ + ... ], 



(44) 
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justify 

= (3.3±0.8) X 10"^ (45) 

as a sensible estimate of the background term in the Onsager coefficient: compare with 
Fig. ^1 The horizontal arrows marked on the right side of Fig. indicate the central 
theoretical value for the amplitude Q. It is evident that the agreement is surprisingly 
good. Indeed, had one been asked to estimate Q from these plots one might have proposed 
Q = (2.7 ± 0.4)xlO~^, again surprisingly close to the theoretical value. Further details of 
this finite-size scaling analysis, including a fit for VVl(T), are presented in Ref. 58. 

Thus we conclude that our simulation data are, in fact, fully consistent with the pre- 
dictions of the theory including the value 0.567 for the exponent ux, and, hence, the result 
xd — 1.0679 for the interdiffusion coefficient itself. It cannot be emphasized too strongly, 
however, that our discussion demonstrates that in the analysis of simulations near critical 
points one needs to account properly for the inevitable finite-size effects and, when theory 
indicates, also for appropriate background contributions typically arising from short-range 
fluctuations.— 

VI. SUMMARY 

We have studied the static and dynamic properties of a symmetric truncated Lennard- 
Jones binary fluid model with ctaa = o"bb = ctab = cr, ^aa = £bb = S^ab = £ and 
masses ttia = = m. This model has a liquid-liquid miscibility gap. We have used a 
combination of semi-grandcanonical Monte Carlo (SGMC) and microcanonical molecular 
dynamics simulations to study both the static and dynamic properties near the demixing 
critical point. The symmetry of the model sets the critical composition at xa = = 1/2. 
We have studied the system at the comparatively high liquid density p* = pa^ = 1, in 
which region the gas-liquid and liquid-solid transitions are far from the temperature range 
of interest. 

The critical temperature has been determined quite accurately as T* = k^Tc/e = 
1.4230 ± 0.0005 using a variety of techniques. Because of the short-range nature of the 
interactions one anticipates that demixing criticality in the model belongs to the three- 
dimensional Ising universality class. All our data for the static properties near the critical 
point strongly support that presumption. 
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We have also presented the first comprehensive study of the dynamic properties of a 
binary fiuid near the critical point. We find evidence for a very weak divergence of the 
shear viscosity, f]{T), near the critical point in accord with expectations. The self-diffusion 
constant D(T) remains finite at the critical point which is consistent with some earlier 
studies. We also find that the Stokes-Einstein relation remains a fairly good approximation 
even within 0.5% of Tc. 

In contrast to the self-diffusion constant, the interdiffusion constant Dab{T) vanishes 
rapidly when T T^. Our analysis of the simulation data supports the various theoretical 
predictions for the critical exponents of all these quantities including the dynamic exponent 
relatio n^ xd = But, even with an accurate knowledge of and of the correlation 

length and concentration susceptibility, it proves essential to consider the finite-size effects 
and allow for background contributions arising from short-range fluctuations, in order to 
properly analyze the data for the interdiffusion coefficient. 

Finally, however, we have not discussed the bulk viscosity, 775 (T), which is expected to 
diverge much more rapidly than the shear viscosity.— That remains a signiflcant task for 
future work. 
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FIG. 1: Probability distributions P{x\) of the relative concentration xa = Na/N of A particles 
for N = 6400 and chemical potential difference A/x = at several temperatures (a) below Tc and 
(b) above Tc, respectively. For clarity many independent data points have been omitted. 
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FIG. 2: Coexistence curve of the symmetrical (truncated) Lennard-Jones binary fluid in the plane 
of temperature T and concentration x\ = Np^/N, for overall density p* = 1.0, the precise choice 
of potentials being given in Eqs. Open circles are the simulation results for a system of 

N = 6400 particles, while the broken curve is only a guide to the eye. The solid curve indicates a 
fit to Eq. ()13() which yields T* = 1.423 as highlighted by the horizontal dot-dashed line. 
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FIG. 3: Finite-size scaling plots of the susceptibility x* for temperatures above using the trial 
values of T* marked in the figure. The Ising values 7 = 1.239, v = 0.629, have been accepted and 
simulation results for x* temperatures T* = 1.45,1.46,1.48,1.50,1.52, and 1.55, are presented. 
Particle numbers from N = 400 to = 12800 are included, as indicated (while the linear dimen- 
sions of the simulation box are L = N^/^a). The dashed lines are guides to the eye: in light of the 
degree of data collapse and the expected scaling function behavior stated in Eq. (|17|) . the estimates 
T* = 1.423 and 1.421 are quite acceptable. 
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FIG. 4: The fourth-order cumulant Ul{T) plotted vs. T for several system sizes, as indicated in 
the figure. The broken horizontal line indicates the value of the Ul at Tc for Ising type systems. 
The vertical line at T* = 1.423 represents our preferred estimate of T*. The smooth curves in the 
enlarged plot (b) are fits to tanh functions. 
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FIG. 5: Plot of the structure factors (a) Snn{Q), (b) Scc{q), for various temperatures, versus 
momentum q. The various curves are shifted up by 0.2 relative to one another for clarity. All 
data refer to a system of = 6400 particles. Inset in part (b) represents an Ornstein-Zernike plot 
which yields estimates for ^ (T) via Eq. . 
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FIG. 6: Plots of (a) the reduced susceptibility x* and (b) the correlation length ^ versus e. Part 
(c) shows the variation of x with ^. The lines represent fits using the anticipated Ising exponents. 
All the data refer to systems of N = 6400 particles. 
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FIG. 7: (a) Log-log plot of the mean square displacements of all the particles versus time with 
to = (ma^/e)"^/^, for systems containing N = 6400 particles, at the critical concentration and the 
seven temperatures indicated. The plots for different T are displaced by factors of 2. (b) Variation 
of the reduced self-diffusion constant D* with temperature. 
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FIG. 8: A log-log plot of the reduced shear viscosity 77* vs. temperature. The line represents a 
least squares fit to the theoretical form with = 0.068 and v = 0.629, yielding an amplitude 
r/o = 3.87 ± 0.3. 
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FIG. 9: Plot of the Stokes-Einstein diameter, d, as defined in Eq. ()29|). vs. temperature. The 
dashed hne serves as a guide to the eye. 
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FIG. 10: Plot of the interdifFusion coefficient ^^^(t) vs. time at tliree different temperatures for 
systems of iV = 6400 particles. The knees visible at short time are due to discrete integration time 
step Ai*. 
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FIG. 11: Log-log plot of the interdifFusion coefficient -D^g as calculated vs. T. The line is a fit to 
the power law I?ab ~ ^^"^^ which yields x^s — 1.6. The data correspond to = 6400. 



33 



2.5 



-ft 



1.5 - 



1.0 - 



0.5 - 



1 — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 



10 C(T} 



2.01% 



5 



5 



J I I L 



J I I L 



J I I L 



J I I L 



J I I L 







0.2 



0.4 



0.6 



0.8 



1 



€ =(T-T )/T 

^ c ^ c 



FIG. 12: Plot of the reduced Onsager coefficient >C(T) vs. T for a system of = 6400 particles. 
Note the "background" contribution and the sharp rise as is approached. The four highest data 
points span the range from 1.9% to 4% above Tc; but the experimentsii probe the range e = 10~^ 
to lO"! 
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FIG. 13: Finite-size scaling plots for the interdiffusional Onsager coefficient jC-l{T) with e = 
(T — Tc)/Tc, y = L/^{T), and trial values for the effective background contribution C^^^ ■ The 
approximate Ising value iy\ = 0.567 has been adopted and, for convenience, we have set = 7 \n 
the abscissa variable, y/{yo + y), that approaches unity when L ^ co. The filled symbols represent 
data at e ~ 4.0 x 10~^ for different system sizes oi N = 400 to 6400 particles and fixed density 
pa^ = 1. The solid arrows on the right hand axis indicate the central theoretical estimate for the 
critical amplitude Q: see text. 
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